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Abstract 



We have performed numerical simulation of a 3-dimensional elastic 
medium, with scalar displacements, subject to quenched disorder. In the ab- 
sence of topological defects this system is equivalent to a (3 + l)-dimensional 
interface subject to a periodic pinning potential. We have applied an efficient 
combinatorial optimization algorithm to generate exact ground states for this 
interface representation. Our results indicate that this Bragg glass is charac- 
terized by power law divergences in the structure factor S{k) ~ Ak^^. We 
have found numerically consistent values of the coefficient A for two lattice 
discretizations of the medium, supporting universality for A in the isotropic 
systems considered here. We also examine the response of the ground state 
to the change in boundary conditions that corresponds to introducing a sin- 
gle dislocation loop encircling the system. The rearrangement of the ground 
state caused by this change is equivalent to the domain wall of elastic de- 
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formations which span the dislocation loop. Our results indicate that these 
domain walls are highly convoluted, with a fractal dimension df = 2.60(5). 
We also discuss the implications of the domain wall energetics for the stability 
of the Bragg glass phase. Elastic excitations similar to these domain walls 
arise when the pinning potential is slightly perturbed. As in other disordered 
systems, perturbations of relative strength 6 introduce a new length scale 
L* ~ beyond which the perturbed ground state becomes uncorrelated 

with the reference (unperturbed) ground state. We have performed scaling 
analysis of the response of the ground state to the perturbations and obtain 
C = 0.385(40). This value is consistent with the scaling relation ( = df/2 — 0, 
where 6 characterizes the scaling of the energy fluctuations of low energy 
excitations. 
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Observation of glassy behavior in flux line arrays in high Tc superconductors [|l| calls for 
a thorough theoretical description of such behavior. In this system the collective pinning 
of the flux line array, rather than the interactions of a single flux line with the disorder, 
can dominate the physics 0. For weak pinning, where dislocations are believed to be 
unimportant at large length scales 0,^, the entire flux line array can be modeled as a 
single medium subject to a pinning potential. Analytic calculations carried out using the 
approximation of linear elasticity and including the effects of the short range order in this 
system indicate that quasi- long range order exists in 3 dimensions The elastic medium 

assumption was justified a posteriori and is further supported by an approximate domain 
wall renormalization calculation p. The structure factor of a topologically ordered system 
was predicted by these calculations to have power law divergences of the form S{k) ~ k^^. 
We will consider only the case of scalar displacements, which also models a charge density 
wave pinned by charge impurities [^. 

We have numerically generated ground states for an elastic medium subject to quenched 
point disorder in the topologically ordered phase. Our results for the coefficient of the 
divergence of S{k) lie between the renormalization group and Gaussian variational method 
results obtained by Giamarchi and Le Doussal [Q. In addition to supporting their analysis, 
we are able to examine the response of the system to changes in boundary conditions and 
pinning potential. By a suitable choice of the boundary conditions we can simulate the 
domain wall of elastic deformations induced by a dislocation loop 0. The energy of the 
domain wall dominates the random part of the energy cost of introducing a single topological 
defect 1^. Our results on the energetics of the domain walls thereby indirectly support the 
analysis carried out by Fisher [Q, which indicated that the this system is marginally stable 
with respect to the introduction of dislocations. The numerically generated domain walls 
were found to have a fractal dimension df = 2.60(5). At large length scales the ground 
state is highly sensitive to small perturbations in the disorder potential, as in spin glasses 
and other disordered systems p-pT|. Perturbations of relative strength 6 in the disorder 
decorrelate the ground state on length scales L* ~ 5"^/^, with ( = 0.385(40). We are able 
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to relate this response to disorder perturbations to the properties of the domain walls. 

We have generated exact ground states for a discrete model whose energy in the contin- 
uum limit is given by 

H = J d^x^ [VM(f)]2 + V{u{x),x) (1) 

with distortions of the medium represented by u{x), which is assumed to be slowly varying 
over the system. The coefficient c is the elastic constant. The potential felt by the medium 
due to the randomly placed impurities is represented by V{u{x),x). In microscopic descrip- 
tions of an elastic system subject to weak disorder there is a length scale below which 
the elastic energy dominates and the medium is ordered. This short range order manifests 
itself here as correlations in the disorder potential of the form V{u,x) = V{u + a,x), with 
a the intrinsic period of the medium. The period of the potential, a, is the lattice spacing 
in a flux line array or the wavelength of a charge density wave. Although this Hamiltonian 
is insufficient to describe the core of a dislocation loop, it can serve to describe the sheet of 
elastic deformations which span the loop, since the approximation of linear elasticity breaks 
down only in the region near the dislocation core. Comparing the ground state of Eq. (|Tp 
in a system of size L subject to periodic boundary conditions to that for the same disorder 
realization with twisted boundary conditions in one direction (i.e., u{x, y, 0) = u{x, y, L) + a) 
allows for the identification of the domain wall which would be caused by a single dislocation 
loop encirchng the system [Q. 

The elastic Hamiltonian, Eq. (|lD, describes a (3 + l)-dimensional interface subject to a 
disorder potential. In this picture, the displacement variable u{x) maps to the height of 
the directed interface. This interface model has a natural discrete representation in which 
the configuration of the interface is specified by the set of bonds it cuts in a 4-dimensional 
lattice ||12[. The bonds of this lattice are assigned weights which directly correspond to the 



disorder potential. The sum of the weights of the bonds which the interface cuts give the 
energy of the configuration corresponding to the disorder energy / d^x V{u,x). In these 
discrete models an effective elastic constant arises from a dependence of the number of 



configurations on the average gradient of tlie interface. Maximal flow algorithms [0, a 
subclass of combinatorial optimization algorithms, allow for the generation of ground states 
of this discrete representation of the interface . 



The lattices numerically studied here are composed of x U nodes, where L is the 
linear size of the elastic medium, and U is the extent of the lattice in the u direction in 



which the displacement variable fluctuates. Unlike the simulation of elastic manifolds [[14 
where the bond weights are non-periodic, there are long range correlations in the disorder 
in the u direction. We generated random integer weights, chosen from a uniform distribu- 
tion over [0, Vmax], independently for each of the forward bonds in a layer of unit cells at 
constant height. The data we will present were obtained with Vmax = 5000, but we have 
verified that our results on the structure of the interface are not significantly altered for 
Vmax as low as 100. Throughout the following discussion, we have normalized the energy 
of the system, so that the effective range of the bond weights is in [0, 1]. This set of bond 
weights on one layer is sufficient to fix the value of the disorder on all of the bonds because 
we require that bonds which differ only by translations along u have the same weight. In 
order to study the universality of the coefficient of the divergence in S{k), we have simu- 
lated interfaces on the simple hypercubic lattice (SHC) as well as the Z-centered hypercubic 
lattice (ZHC) [T^. For the SHC lattice u was chosen to be along the (1111) crystallo- 
graphic axis. The elementary bonds in this lattice are: {x,y,z,u) = ±(0, 1/2, 1/2), 

±(0,-^2/2,1/2,1/2), ±(V2/2, 0,-1/2, 1/2), ±(-^2/2, 0, -1/2, 1/2). In the ZHC lat- 
tice we considered the following 12 bonds extending from each node: [x, y, z, u) = 
(±V2/2,0,0,±V2/2), (0,±V2/2,0,±V2/2), (0, 0, ±^2/2, ±^2/2) These lattices are 
the natural extensions of the two types of lattices used in simulations of the ground state 
of a 2-dimensional elastic medium [|T^ . Both types of lattices were simulated using periodic 
boundary conditions in the transverse directions. In addition, in the ZHC lattice the ground 
state was computed for each realization of disorder with twisted boundary conditions in one 
direction. 

We have developed a custom implementation of the push-relabel maximal fiow algorithm 



18| optimized for application to the regular lattices considered here. Our modifications to 
the Push-Relabel algorithm reduced the memory requirements by nearly a factor of ten, 
allowing for simulation of systems composed of up to approximately 6 x 10^ nodes using 
less than 512MB. The primary modification involves computing nearest neighbor relations 
as needed rather than storing this information. Overhangs in the interface are precluded 
by assigning a large weight to backwards arcs fl^. Since these backwards arcs have an 



effectively infinite weight, they cannot be be part of the minimal cut. Thus our algorithm 
can operate without storing their weight, provided that fiow is always allowed to move along 
the backwards arcs. These modifications increased the running time of the algorithm, but 
obtaining the ground state for each realization of disorder still took less than one hour of 
processor time on a single 400 MHz Pentium II CPU for the largest system sizes studied. 
The memory requirement is linear in the number of nodes A^; the processor time was found 
to scale approximately as A^^'^, compared with the worst case bound of A^^ |T^. 

The modest computational requirements of this algorithm have allowed us to average 
the properties of the ground states for a variety of system sizes over a large number of 
disorder realizations. In addition to generating the value of the minimal energy, the algo- 
rithm produces the configuration of the interface. The interface can be then represented by 
u{x), which is defined on the 3-dimensional lattice formed by projecting the interface along 



the u direction |jT9[. Due to the periodicity of the disorder, the energy is invariant under 
global translations of integer multiples of a in the displacement variable u{x) — u{x) + na. 
Considering the set of the forward bonds cut at each location provides a characterization 
of the configuration equivalent to measuring the gradient of the interface. This represen- 
tation of the interface is useful when comparing different ground states since it insensitive 
to global shifts of u. For the SHC lattice, we have generated at least 10^ realizations of 
disorder for systems of size L = 8, 16, 24, 32, 40, 48, 60, 80. We chose the extent of the lattice 
in the displacement direction u to ensure that the boundaries of the system do not affect 
the ground state. For the SHC lattice U = 20 was sufficient. Our simulations for the ZHC 
lattice were more extensive, with at least 10^ realizations for systems, subject to both peri- 



odic and twisted boundary conditions, of size L = 8, 16, 32, 48, and at least 10^ realizations 
for systems of size L = 64,80. The largest systems here required U = 12 to prevent the 
configuration from being affected by the boundaries of the lattice in the u direction. 

We have examined the displacement correlations of the minimal energy configurations by 
computing the disorder averaged structure factor S{k) of the displacement variables. This 
allows us to more clearly distinguish the large length scale behavior; direct measurement of 
the width is more difficult to analyze due to finite size effects. The orientationally averaged 
structure factor S{k) has been obtained by averaging the value of S{k) over radial bins 
of size Ak = 0.025, and is presented in Fig. |l|. The error bars represent the fluctuations 
of S{k) within each spherical shell. These fluctuations, which measure the anisotropy of 
the structure factor, generally decrease with decreasing k; for k < 0.75, these fluctuations 
saturate at a value comparable to the statistical fluctuations in S{k) indicating the range in 
k where the system is isotropic within the statistical fluctuations. To extract the coefficient 
of the leading order divergence we have fit k^S{k) with the form 

P'S{k) =A + Bk (2) 

over the region k < 0.5. The leading order term of S{k) ~ k~^ indicates the quasi-long 
range order of the ground state. For the ZHC lattice we obtain k^S{k) = 1.08(5) + 1.02(6)A;; 
for the SHC lattice we obtain k^S{k) = 1.01(4) + 0.46(4)A;. The error estimates on the 
parameters in these fits represent the statistical uncertainty over the given fit range and 
systematic errors arising from the choice of the cutoff in k. Giamarchi and Le Doussal have 
applied a renormalization group technique to order e = 4 — d to this system and obtained 
A = 1.0 0. They have also carried out a Gaussian variational calculation and obtained 
A = 1.1. In both of these approximations the value of the coefficient is universal. 

The form Bk~^ of the leading order corrections to S{k) can be obtained by a renormal- 
ization group calculation. Upon renormalization, the periodic pinning potential V{u{x),x) 
introduces a new term into the Hamiltonian of the form fl{x) •V'u(x) as when d = 2 P,pD|. 
This random tilting field is short range correlated with /2(x) -fliy) = g'^6^{x — y). Unlike the 
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case where d = 2, the strength of this field g"^ undergoes only a finite renormalization for 
d = 3 [^]. The effects of this term can be determined from the effective small length scale 
Hamiltonian 

Heff = j d^x ^[\/u{x)f + •VM(f) , (3) 

which ignores the periodic pinning potential. Solving the equations of motion and averaging 
over realizations of the tilting field predicts corrections of the form Bk~'^ with B = g'^/c^ at 
T = 0. In order to consider the effects of the pinning potential and the tilting field separately 
we relied on the fact that the renormalization group fiow of V{u{x),x) is unaffected by the 
presence of a non-zero g PJ20||. The results of this analysis can be confirmed by examining 



the stability of the functional renormalization group fixed point obtained by Giamarchi and 
Le Doussal in an e = 4 — expansion [Q. The most slowly decaying perturbation to the 
fixed point decays ~ L"*^, implying corrections to the structure factor of the form ~ 
This order e calculation also predicts the form of the corrections which are observed in our 
simulation data. 

Naturally, real-space measurements of the width must be consistent with the structure 
factor. We measured the disorder averaged squared width 

ty^ = <m2> - <M>2 , (4) 

where < > denotes a spatial average over the system. These data, shown in Fig. |^, have been 
fit using the real space version of Eq. (2): w"^ = a + bln{L) + c/L . The constant term arises 
from the short wavelength fiuctuations, while the second and third terms arise from the 
and terms in the structure factor respectively. The real space coefficient b is related to 
the leading order behavior of S{k) hjb = A/A-n-"^. This three parameter fit gives estimates for 
the values of A and B, the coefficients describing the long wavelength form of the structure 
factor, which are consistent with that obtained by fitting S{k) at a single system size. If 
the form of S{k) depended on L, then this consistency would not be maintained. Direct 



comparison of the structure factor for various system sizes also demonstrates that S{k) 
exhibits negligible system size effects, except for the change in k^in = L/2tt. 
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Other measures of the displacements provide us with additional information on the 
structure of the ground state. The disorder averaged extremal displacement difference, 
AiJ = Umax — Umin , Fig- 0, was fouud to grow logarithmically with system size for both 
lattice types. We computed least squares fits of the form AH = d + b ln{L) to obtain 
bzHC = 0.76(1) and hsuc = 0.70(1). The coefficients of the logarithmic term differ by less 
than 10% for the two lattices studied here, suggesting that this measure of the system is 
weakly, if at all, dependent on the lattice discretization of the medium. This logarithmic 
growth is consistent with the following picture of the ground state structure developed by 
Fisher |12||. At each length scale R = b, b^, b^, . . . the displacement undergoes one shift of 



amount ±a. Furthermore the sign of the displacement shift is random at each scale, leading 
to the logarithmic growth of the squared width. When traversing from the minimum to the 
maximum the signs of the displacements are strongly correlated, leading to a coherent sum, 
and a logarithmic dependence on the system size for the extremal differences results. 

We have also determined the effect of coarse-graining the displacement variable. The 
coarse grained displacement is defined as the average of u, 

UR{y) = 4? / u{x) , (5) 

over flR^y), a cube of size R centered at the point y. We measured the fluctuations in these 
coarse grained height variables | Am^jP = {uR{y) — UR{y + b)Y, with b the vector between the 
centers of the cubes which touch at one corner. This spatial averaging procedure is similar to 
a real space renormalization transformation. Villain and Fernandez have explicitly carried 
out a real space renormalization calculation for a 3-dimensional elastic medium with cubic 
symmetry [^. Their calculations indicate that IAmr^ has a finite limit as R,L oo. 
We have directly measured IAm^jP for the ZHC lattice (Fig. |^). The coarse grained height 
fluctuations are related to the structure factor by: 

with G{k) = J^^d^x e*'^'^. In the infinite volume limit IAmrP depends only on the leading 
order behavior, Ak~^, of the structure factor and the limit of the ratio R/L. We have 



numerically evaluated the right hand side of Eq. (j^) in this limit to obtain the infinite size 
limit presented in Fig. ^. For finite sized systems the sub- leading order corrections to S(k) 
contribute to the coarse grained height fiuctuations. These terms lead to the divergence of 
lAu/jp as R/L 0. The dominant corrections arise from the Bk'"^ corrections to S{k) seen 
in the structure factor data, but decay as L~^. The data are consistent with convergence to 
a finite limit as L ^ oo. 

For the ZHC lattice, we have also investigated the behavior of the system subject to the 
twisted boundary conditions defined previously. For each realization of disorder we have 
compared the ground state energy with periodic boundary conditions, Ep, to that obtained 
with twisted boundary conditions along one of the lattice directions, Ef. This allows us to 
investigate the properties of the excitations induced by the change in boundary conditions. 
We identify the energy difference -E_dvk = Ef — Ep with the energy of the domain wall. 
The domain wall is identified by the set of bonds that the interface intersects in one set of 
boundary conditions but not the other. Even though the domain wall could be identified by 
examining the values of Up{x) and Ut{x), it is more efficient to identify the domain wall by 
examining the sets of cut bonds for each boundary conditions. The cut bonds are projected 
along the u direction for each boundary condition. The domain wall is then the symmetric 
difference between these two sets of bonds. Without directly simulating a dislocation loop 
itself, we are able to investigate the properties of the domain wall induced by the introduction 
of a single loop encircling the system. 

The energetics of these domain walls dominate the random part of the energy cost of 
introducing a dislocation loop into an elastic medium P]. The mean energy difference Edw 
grows linearly with the (linear) size L of the domain wall (Fig. ^), a result consistent with 
the scaling of the elastic contribution to the energy and the statistical symmetry of the 
disorder potential of the continuum model We have also analyzed the variance a'^^ED]^) 
of the distribution of domain wall energies (Fig. No single power law fit of the form 
a'^i^Eow) oc can adequately fit our data over the range of sizes simulated. However the 
data are well described by the empirically determined form a'^{EDw) = O.OSIL^ + 0.24L 
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displayed as the solid line in Fig. ^. The finite size correction leads to a size dependent 
effective value of the exponent, Oeff, characterizing the scaling of the sample to sample 
fiuctuations in the energy of low energy excitations. Depending on the lower limit imposed 
on the fit, single power-law fits a'^^Ecw) ~ L"^^ give 0.85 < 9eff < 0.92 over this range of 
sizes. The distribution of domain wall energies depends on and a{ED]y) only through 
the combination e = {Edw — Enw) / o-^Edw), as can be seen in Fig. 0. This collapsed 
distribution has more highly weighted tails than a unit normal distribution. The frequency 
with which negative energy domain walls were observed in the simulations is significantly 
higher that what one expects for a Gaussian distribution with the measured mean and 
standard deviation at each system size (Fig. §). This behavior also occurs for high energy 
domain walls; to within the statistical uncertainty the distribution is symmetric about its 
mean value. Our data are consistent with both E^w, and a{Eow) increasing linearly with 
L for large systems. Fisher's argument assumed this behavior of the domain wall energetics 
in his domain wall renormalization calculation indicating the marginal stability of the Bragg 
glass phase 0. 

Because of the balance between the elastic energy scale and the scale of the energy 
fiuctuations due to the disorder potential, the domain walls are highly convoluted and ex- 
pected to have a fractal dimension df between 2 and 3. Similar to the approach used for 



a 2-dimensional elastic medium we have measured the size of the domain wall by 

counting the number of bonds A^^f, in the wall. The data for the area of the wall as a function 
of system size, averaged over disorder, can be fit by a simple power law, Nf, ~ L'^f , with 
the fractal dimension of the domain wall df ~ 2.60, shown in Fig. ^(a). This fit has been 
taken after excluding the smallest system size L = 8. However, the form of the residuals, 
see Fig. ^(b), indicate the presence of sub- leading order corrections, suggesting that this 
may underestimate the value of df. In order to verify this, we have also fit the whole range 
of data using the form A*";, = aL'^f + bL^. The second term arises from the effectively two 
dimensional nature of the domain walls at small length scales. This three parameter fit gives 
df = 2.65; thus we conclude that our systematic errors in estimating the fractal dimension 
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are approximately 0.05. 

In order to investigate other low energy excitations of this system, we have examined the 
sensitivity of the ground state to perturbations in the disorder potential. Similar studies 
have been carried out for disordered systems such as spin glasses 0, (1 + l)-dimensional 
directed polymers in random media [|TU|, and 2-dimensional elastic media |]TT|. In all of these 



disorder dominated systems, perturbations of relative strength S in the disorder potential 
introduce a length scale L* ~ 5"^/^, ( = df/2 — 9, beyond which the ground state becomes 
uncorrelated with the reference ground state. The exponent ( characterizing the sensitivity 
of the ground state is referred to as the chaos exponent. These studies have been done 
by comparing the ground state for two correlated choices of the disorder potential. In our 
simulations, we have obtained the ground states for the two pinning potentials V^{u{x),x) = 
b{u{x),x) ± d{u{x),x), with both terms periodic in the u direction. The constant part 
of the potential b{u{x),x) was an integer chosen uniformly from [1000,2000]. The term 
which generates differences between the realizations, d{u{x),x), was chosen uniformly from 
[—dmax/'^jdmax/'^]- The parameter 6 = dmaxf^^^O characterizes the relative strength of 
the perturbations. This prescription was chosen to ensure that for a fixed value of 6 the 
distribution of the bond weights is the same for both realizations of disorder. Our simulations 
include values of S ranging from 0.01 to 0.75, with at least 500 independent realizations of 
disorder at each 6 and L. By performing scaling analysis of both the energetic and structural 
correlations between the ground states for these realizations of disorder we can extract the 
value of the chaos exponent. 

We have found that both the energetic and structural correlations are governed by the 
same length scale. First we calculated the domain wall energy E^^^ for the two disorder 
realizations V^, which can then be used to compute the domain wall energy correlation 
function ISIl 



The simple scaling form G = f{6L'^) describes our data well (Fig. |ToD. We found reasonable 
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data collapse for ( = 0.38(4), taking into account the statistical errors. The value of the chaos 
exponent can be related to the domain wall fractal dimension and the energy fluctuation 
exponent by a simple scaling argument as in the case of spin glasses |p. The perturbations 
introduce a random change in the energy of the domain wall of order 5L'^f^'^ because the 
perturbations are uncorrelated with the location of the domain wall. The typical fluctuations 
in the domain wall energy scale as L^. When these energy scales become comparable, at 
a length scale L* ~ 5^^^^, the domain wall energies become uncorrelated. When using the 
effective value of the energy fluctuation exponent, 9eff ~ 0.9, this scaling relation holds to 
within 5% accuracy. 

We can understand the structural deformations induced by the bond perturbations by 
reasoning similar to that for domain wall correlations. Here we consider the differences in the 
ground states of the system (with periodic boundary conditions) due to the changes in the 
pinning potential. Again, the ratio of the energy change due to the random perturbations, 
and that of the fluctuations in the energy landscape, SL'^f^'^/L^, determines the behavior 
of the system. In response to Zhang's simulation of directed polymers in random media, 
Feigel'man and Vinokur had argued that the probability of a positional excitation grows 
linearly with SL"^ for small values of the perturbation strength [^. Following their argument. 



we expect the probability of a change in the displacement variable at a single location to grow 
linearly with 6L''. Unlike the case for the directed polymer, the magnitude of the differences 
is bounded for periodic pinning since excitations with |Aii| > a probe the same energy 
landscape as those with |Am| < a. Thus, the spatially and disorder averaged mean squared 



displacement difference x = <(m^(x) — u^{x)y> ~ for small perturbations. For large 
values of 6 the ground states are completely decorrelated, and we expect that x ~ ln{L). 
We have found that the data for x collapse according to the scaling form x = fi^L'^), 



with C = 0.39(2) (Fig. |Tl]). Our data are consistent with the results obtained by the 
scaling argument in both limiting cases. Before computing x, we made the transformation 
— i> + na, where n is the integer which maximizes the number of locations at which 
u^{x)—u~{x) = for each realization. This transformation minimizes x over the discrete set 
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of global translations which leave the energy of the medium invariant. Our scaling ansatz is 
significantly different from that proposed for the 2-dimensional elastic medium x = ^^^^/(-^) 
||1 1|| . This form cannot adequately collapse our data over the range of parameters we have 
simulated. 

Implicit in this discussion is the assumption that both the perturbation induced defor- 
mations, and the boundary condition induced domain walls are characterized by the same 
fractal dimension. Even for relatively small perturbations in the disorder, the deformations 
are typically composed of a set of disconnected clusters. We have directly measured these 
clusters' fractal dimension. The size of a cluster R is defined as the average of the sides of 
the bounding box which encloses the cluster and is measured in units of the lattice spacing. 
Our algorithm identifies the sets of nodes on which u~^{x) — u^{x) ^ after performing the 
translation which minimizes x- The surface area of a cluster s is the number of nodes with 
neighbors not in the cluster. We have collapsed the data for 6 = 0.05 using the finite size 
scaling form s = L'^f f{R/L) (Fig. |T2|). We expect that the scaling function f{R/L) should 
have the form / ~ {R/LYf in the region R/ L « 1, R » 1, but this regime is not clearly 
visible in Fig. |T2| due to lattice and finite size effects. Despite this, the best collapse of the 
data, for 0.5 < R/L < 1.0, provides an estimate df = 2.65(10) which is consistent with 
the estimate for the domain wall fractal dimension. Similar analysis at other values of S 
provide equivalent values for the cluster fractal dimension. The anomalous data at ^ 1 
arise from the rare clusters which span the system in all directions. This scaling also breaks 
down for clusters with R < 8, where lattice effects make the surface effectively 2-dimensional 



(Fig. |T^). The equality of the fractal dimension of the boundary condition induced domain 
walls and bond perturbation induced deformations can be justified by a simple argument. 
For small values of the disorder perturbation parameter, the cluster boundaries lie in regions 
where there is a small energy cost to deforming the medium. If one considers only a small 
volume containing a portion of the cluster boundary, the structural difference is the same 
as would be caused by the change from periodic to twisted boundary conditions on that 
volume. Thus both the deformations induced by small changes in the disorder potential 
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and those caused by a change in boundary conditions should be characterized by the same 
fractal dimension. Despite the fact that the scaling regime is inaccessible due to the limits 
on the size of systems studied, our data are consistent with the conclusion that the frac- 
tal dimension of the clusters which compose the deformations is the same as that of the 
boundary condition induced domain walls. 

We have performed extensive numerical simulation of a model 3-dimensional elastic 
medium subject to quenched disorder with scalar discrete displacements. Our results for 
the structure in the Bragg glass phase indicate that the structure factor has divergences of 
the form S{k) ~ Ak~^. Our results for the coefficient A fall between the approximate values 
A = 1.0, and A = 1.1, obtained via a renormalization group and a replica approach [Q. 
The observed energetics of the boundary condition induced domain walls indirectly support 
arguments for the stability of the Bragg glass phase. These domain walls correspond to 
the elastic deformations due to the introduction of a single dislocation loop winding around 
the system. Our data are consistent with the hypothesis that the mean energy and the 
energy fluctuations of a section of domain wall both scale linearly with the linear size of 
the section for large sizes. This balance is a crucial element of the analysis carried out 
by Fisher indicating the marginal stability of the Bragg glass phase to the introduction of 
dislocations f^. We are also able to measure the spatial structure of these domain walls and 
obtain their fractal dimension df = 2.60(5). We have observed that random changes in the 
disorder potential of relative strength 6 decorrelate the ground state on length scales larger 
than L* ~ S^^^'^ with = 0.385(40). The properties of the domain walls and this sensitivity 
to disorder perturbations can be related to each other by the scaling relation ( = df /2 — 6, 
where 9 characterizes the fluctuations in the low energy excitations. 
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FIGURES 

FIG. 1. Numerically calculated structure factor for a pinned elastic medium, averaged 
over disorder realizations. Data is shown for: (a) the Z-centcrcd hypercubic lattice for L=80 
and 1000 realizations, and (b) the simple hypercubic lattice, L=80, 1000 realizations. The 
data have been coarse grained into bins of size Ak = 0.025, and only every third data point 
is displayed. The error bars represent the fluctuations in the values of S{k) within each 
bin, and hence are a measure of the anisotropy of S{k). For both lattices a fit of the form 
k^S{k) = A + Bk has been taken over the region k < 0.5. The resulting least squares fits 
are shown as solid lines. For wave vectors in this region the anisotropy of the structure factor 
is comparable to the size of the statistical fluctuations in the value of S{k). The coefficient of 
the leading order divergent term, Ak~^, can be extracted from these fits: Azhc = 1.08(5), 
AsHC = 1-01(4), where the errors include both the statistical errors and our estimate of the 
systematic errors. 

FIG. 2. Disorder averaged roughness of the interface representation of the elastic medium. 
The statistical errors at each size are comparable to, or smaller than the plot points. Fits of the 
form a + b ln{L) + c/L, which include the form of the finite size corrections indicated by the 
structure factor data, are represented by the solid lines. The coefficient h in these fits is related to 
the coefficient of the divergence of the structure factor by 6 = A/Att^. These fits provide values of 
AsHC = 1-03(6), Azhc = 1-02(1). 

FIG. 3. Disorder averaged extremal displacement difference, AH = Umax — Umin , as a function 
of system size. The statistical errors at each size are smaller than the plot points. Over this 
range of sizes the behavior is logarithmic, as indicated by the solid lines which are fits of the 
form AH = a + b ln{L). The least squares values for the coefficient of the logarithmic term are 
bzHC = 0.76(1), bsHC = 0.70(1). 
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FIG. 4. Behavior of the coarse grained displacement differences. | An| is the difference between 
the average displacement of the medium between two adjacent cubic regions of size R, averaged 
over disorder. The statistical errors at each size are smaller than the plot points; the dashed lines 
serve as guides to the eye. The solid line represents the L ^ oo extrapolation for the coarse grained 
displacement differences and was calculated using the long wavelength behavior of the structure 
factor. 

FIG. 5. Average energy difference between the ground state of the system with periodic and 
twisted boundary conditions. The data are well fit by Edw = 0.395(1)L. 

FIG. 6. Variance of the energy difference Edw between the ground states with periodic 
and twisted boundary conditions. The fit, cr^ = 0.031(1)L^ + 0.24(1)L, includes an empirically 
estimated correction to scaling. 

FIG. 7. Data collapse for the distribution of e = {Edw — Edw) /(^ {Edw) ^ the normalized 
energy difference between the ground states of the system with twisted and periodic boundary 
conditions. The points represent the observed frequencies of e over bins of size 0.2 for at least 10^ 
realizations at each system size. 

FIG. 8. Frequency with which negative energy domain walls were observed. The solid line 

represents the expected probability given the assumption that the distribution is a Gaussian with 
the observed mean and standard deviation. 

FIG. 9. Structural difference between systems with periodic and twisted boundary conditions. 
Part (a) displays the area Ni, of the domain wall as a function of the system size L. The statistical 
errors at each size are smaller than the plot points. The solid line represents a fit over the region 
L > 10 with the form /(L) = aL'^f . The least squares fit provides an estimate of the domain wall 
fractal dimension df = 2.60. However the form of the residuals for this fit, part (b), suggest the 
presence of small sub-dominant corrections. 
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FIG. 10. Boundary condition induced domain wall energy correlation function G. The scaling 
variable is a combination of the system size L and bond perturbation strength 5. The typical 
magnitude of the of these perturbations ranges from 5 = 0.01 io 5 = 0.75 when measured in units 
of the range of the (unperturbed) pinning potential. Each data point represents an average over at 
least 500 realizations of disorder. The data collapse provides an estimate for the chaos exponent 
C, = 0.38(4); the error was estimated by determining the range in where the collapse is reasonable. 

FIG. 11. Configurational differences induced by perturbing the disorder potential randomly 
at each location with relative strength 5. The response function, % = '^xi'^^'^i^) ~ is 
the spatially averaged squared displacement difference between the reference and perturbed ground 
states. The data have the scaling form x = f{^L,^)^ with C, = 0.39(2). 

FIG. 12. Scaling of the surface area of the singly connected components of deformations 
induced by perturbations of relative strength 5 = 0.05 in the disorder potential. The data collapse 
for df = 2.65, with an uncertainty of ±0.1. As a guide to the eye, the df = 2.65 power law behavior 
of the average surface function of the cluster's radius, R, is indicated by the solid line. 

Because each realization of disorder generates numerous clusters of various sizes, the statistical 
errors for the average surface area are smaller than the plot points. 

FIG. 13. Relationship between the surface area and size of small connected clusters of deforma- 
tion induced by small changes in the disorder potential. For R < 8 their surfaces are approximately 
two dimensional and independent of the system size, as indicated by the solid line representing 
s = R^. This behavior does not persist to larger clusters. In order to illustrate the cross over, the 
large cluster behavior s ~ R^-^^ is represented by the dotted line. The statistical errors for each 
data point are smaller than the plot points. 
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Fig. 1: D. McNamara, A. A. Middleton, C. Zeng "Simulation of the Zero ..." 





Fig. 4: D. McNamara, A. A. Middleton, C. Zeng "Simulation of the Zero ..." 
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Fig. 8: D. McNamara, A. A. Middleton, C. Zeng "Simulation of the Zero ..." 
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Fig. 11: D. McNamara, A. A. Middleton, C. Zeng "Simulation of the Zero ..." 




Fig. 12: D. McNamara, A. A. Middleton, C. Zeng "Simulation of the Zero ..." 




